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We investigate the application of the Density Matrix Renormalization Group (DMRG) to the 
Hubbard model in momentum-space. We treat the one-dimensional models with dispersion relations 
corresponding to nearest-neighbor hopping and 1/r hopping and the two-dimensional model with 
isotropic nearest-neighbor hopping. By comparing with the exact solutions for both one-dimensional 
models and with exact diagonalization in two dimensions, we first investigate the convergence of the 
ground-state energy. We find variational convergence of the energy with the number of states kept 
for all models and parameter sets. In contrast to the real-space algorithm, the accuracy becomes 
rapidly worse with increasing interaction and is not significantly better at half filling. We compare 
the results for different dispersion relations at fixed interaction strength over bandwidth and find 
that extending the range of the hopping in one dimension has little effect, but that changing the 
dimensionality from one to two leads to lower accuracy at weak to moderate interaction strength. In 
the one-dimensional models at half-filling, we also investigate the behavior of the single-particle gap, 
the dispersion of spinon excitations, and the momentum distribution function. For the single-particle 
gap, we find that proper extrapolation in the number of states kept is important. For the spinon 
dispersion, we find that good agreement with the exact forms can be achieved at weak coupling 
if the large momentum-dependent finite-size effects are taken into account for nearest-neighbor 
hopping. For the momentum distribution, we compare with various weak-coupling and strong- 
coupling approximations and discuss the importance of finite-size effects as well as the accuracy of 
the DMRG. 

PACS numbers: 71.10.Fd, 71.27.+a 
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I. INTRODUCTION 

Many renormalization schemes are carried out in mo- 
mentum space and involve integrating out degrees of 
freedom using a momentum cutoff. For example, Wil- 
son's numerical Renormalization Group (RG)u imple- 
ments this program using a mapping of momentum shells 
to an effective lattice model. The renormalization process 
is carried out by successive numerical diagonalization of 
a finite system and energetic truncation of the Hilbert 
space. While this lattice model corresponds to succes- 
sive momenta or, equivalently, energy scales, its form is 
similar to that of a strongly correlated lattice model. 

Attempts at applying a real-space version of the Wil- 
son procedure to short-range quantum lattice models 
such as the Heisenberg or the Hubbard model were not 
successful, however, because successive lattice points do 
not correspond to different energy scales. |-.jrhe Density 
Matrix Renormalization Group (DMRG)cH3 overcomes 
these limitations by carrying out the renormalization on 
a subsystem. The truncated basis is formed by pro- 
jecting the state of the entire system onto the subsys- 
tem using the reduced density matrix rather than se- 
lecting states energetically. This method has been very 
successful at treating low-dimensional quantum lattice 
models with open boundary conditions and short-range 
couplings. However, for longer-range off-diagonal inter- 
actions, higher dimensional systems or lattices with peri- 
odic boundary conditions, this real-space formulation of 



the DMRG is much less successful. In addition, it loses 
sight of an energy or momentum-based classification of 
the relevant degrees of freedom. 

A potential way of overcoming this limitation for itin- 
erant electron systems is to apply the DMRG ideas to 
the momentum-space formulation of the Hamiltonian. 
This approach has a number of potential advantages 
over the real-space approach. First, since the single- 
particle basis in momentum space is explicitly transla- 
tionally invariant, momentum is a conserved quantum 
number. Use of this momentum quantum number re- 
duces the size of the Hilbert space in the diagonaliza- 
tion. Second, momentum-dependent quantities such as 
the momentum distribution or the dispersion of excita- 
tions can be directly calculated. Third, the kinetic en- 
ergy term is diagonal so that varying the dispersion by, 
for example, changing the range of the hopping, is easy 
to do. 

Attempts to formulate a numerical renormalization 
group procedure for quantum lattice systems in momen- 
tum spacaj predating the DMRG were not particularly 
successful - this was one of White's motivations for turn- 
ing to real space and formulating the DMRG. Shortly af- 
ter the development of the DMRG in real-space, White 
attempted to use DMRG methods on the momentum- 
space formulation of the Hubbard model. He calculated 
the ground-state energy in one and two dimensions at 
intermediate couplings, but found that the energies ob- 
tained were not significantly better than those obtained 



2 



by other variational methodsB 

Independently, Xiang developed a similar technique 
and applied, it to the Hubbard model in one and two 
dimensions.B In this work, Xiang outlined an efficient im- 
plementation of the DMRG in momentum space. He de- 
veloped a factorization of the Hubbard interaction that 
reduces the number of term from N^, where N is the 
number of single-particle Bloch wavefunctions in the lat- 
tice, to 6N. He also pointed out some features of the 
algorithm that need to be carefully considered in mo- 
mentum space: Since the interaction is highly non-local, 
there is no natural ordering of the single-particle states; 
the choice of the ordering can, however, have an effect 
on the performance of the DMRG algorithm. In addi- 
tion, there is no well-defined infinite-system algorithm, 
so that care must be taken in how the lattice is built 
up initially. Care must also be taken in this initializa- 
tion procedure so that states with a sufficient spread in 
momentum quantum numbers are kept. One possible 
outcome of an inadequate initialization procedure is con- 
vergence to a state other than the true ground state. 

Xiang investigated the performance of the algorithm 
for various interaction strengths for the one-dimensional 
Hubbard model with 16 sites at half filling, and for the 
two-dimensional model on system sizes ranging from 4x4 
to 12 X 12 for various band fillings. He found that the 
convergence of the method depends strongly on the in- 
teraction strength, U. The method is exact for U = 
since the Hamiltonian is diagonal and the convergence 
becomes rapidly worse with increasing U. In one dimen- 
sion, he compared with real-space DMRG calculations 
and found that the error in the ground-state energy was 
higher than the real-space calculation for both weak and 
intermediate interaction strengths {U/t — 1 and 4, with t 
the hopping matrix element), with 5% error for U/t — 4. 
In two dimensions, he compared with exact diagonaliza- 
tion for a 4 X 4 system, cluster diagonalization on a 6 x 6 
system, and quantum Monte Carlo and stochastic diag- 
onalization calculations on 4 x 4, 6 x 6 and 8x8 sys- 
tems. The relative errors increased rapidly with U for 
the 4x4 system. The variational energies were com- 
parable to those of the stochastic diagonalization and 
quantum Monte Carlo methods (for which the energy is 
non-variational) for larger system sizes. The variational 
bounds for the energy were slightly higher than stochas- 
tic diagonalization for the 6x6 lattice and slightly lower 
for the 8x8 lattice. In comparing the performance in one 
and two dimensions, Xiang pointed out that the accuracy 
for given number of states kept, U/t, and band filling for 
16 site systems was better in two dimensions than in 
one, leading him to speculate that the momentum-space 
method becomes more accurate as the dimensionality is 
increased. 

Our purpose in this work is to explore more fully both 
the convergence properties and the application of the 
momentum-space formulation of the DMRG to the Hub- 
bard model. In one dimension, we take advantage of the 
existence of exact solutions for two choices of the disper- 



sion, corresponding to nearest-neighbor hopping and 1/r 
hopping, to systematically investigate the dependence of 
the convergence of the ground-state energy on interac- 
tion strength, band filling, and number of density-matrix 
eigenstates kept, m. We investigate the regularity of 
the convergence with m and discuss schemes to extrap- 
olate in TO in order to obtain more accurate energies. 
We reexamine the relative convergence for the one and 
two-dimensional models with a view to understanding 
the utility of the momentum-space DMRG for higher- 
dimensional systems. 

While the ground-state energy is useful for determin- 
ing variational convergence, it does not directly provide 
much useful information about the physical behavior of 
the system. We therefore investigate some physically use- 
ful quantities, the quasiparticle gap, the momentum dis- 
tribution and the dispersion of spinon excitations, for the 
one-dimensional models and compare to exactly known 
results and perturbative approximations, where appro- 
priate. Our calculations of the momentum distribution 
and the spinon dispersion for the 1/r-hopping model are, 
to our knowledge, the first independent numerical calcu- 
lations of these quantities. 

The layout of the remainder of this paper is as fol- 
lows: In section II we discuss the model systems and 
their basic properties. Our DMRG method is described 
The convergence and accuracy of the 



in section III 



momentum-space DMRG are discussed in section IV 



We study the dispersion of spinon excitations and the 
momentum density distribution of one-dimensional Hub- 
bard models in sections ^ and VI , respectively. We dis- 
cuss the prospects for momentum-space DMRG in the 
final section. 



II. MODEL 

The Hubbard modelB is defined in a general form by 
the Hamiltonian 



(1) 



where (cia) creates (annihilates) an electron with spin 
a in the Wannier state on lattice site i with position r^, 
riiCT — cl^Ci^ denotes the particle number operator on 



' rj) is the transfer integral between site 



site i, tij — tiri - 
i and j, and U is the energy cost due to the Coulomb 
repulsion of two electrons on the same site. In this paper, 
all energies are measured in units of t — 1. 

Using the relation (Fourier transformation) between 
Wannier states at site i and Bloch states with momen- 
tum k 
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E 



(2) 



where N is the number of lattice sites, the Hubbard 
Hamiltonian with translationally invariant transfer inte- 
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grals is written in momentum space as 

Cp_qT4+qiCfciCpt , (3) 



N ^ 

p,k,q 



where Uka = c^^c^a and 



(4) 



is the energy dispersion of the electrons. 

The kinetic energy of Eq. (^) consists only of diagonal 
terms with dispersion e(fc), so that the momentum-space 
DMRG method is trivially exact for U — 0. Moreover, 
it can be easily applied to different non-interacting dis- 
persions e{k) corresponding to different lattice geome- 
tries and hopping ranges. In this paper, we apply the 
momentum-space DMRG to the following three differ- 
ent models (here and in what follows, we take the lattice 
constant to be unity and N to be even): 

(i) The one-dimensional Hubbard chain with nearest- 
neighbor hopping amplitude tj+ij = —te^'^. The disper- 
sion relation is given by 



e(fc) = -2icos(fc - (/)) 



(5) 



with k = 2nn/N and n = -7V/2 + 1, • • •,7V/2. The band- 
width is W = At. Here a flux N(t), measured in units of 
the flux quantum (j)o = hc/e and equivalent to a twiste' 
boundary condition, is threaded through the systemS 
It enables us to calculate the ground-state energy as a 
function of flux (/). We will use nonzero values of 4> later 
to interpolate momentum-dependent quantities (e.g. the 
momentum distribution function) to arbitrary values of 
the momentum on a finite system. 

(ii) The one-dimensional Hubbard chain with long- 
range hopping amplitude 



with 



tlm - {-'^ij-TTi T 

d[L — to) 
N . \T:{l-m) 



d{l — to) = — sin 

TT 



N 



(6) 
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Since d{l — to) is antisymmetric under the permutation 
of I and TO, the hopping matrix element has to be purely 



imaginary to guarantee that tin 



t 



In the thermody- 



namic limit {N — > c»), the hopping decays proportionally 
to the inverse of the distance r = |Z — to| ("1/r-Hubbard 
model"). The dispersion relation is linear and is given 
by 



s{k) = tk 



(8) 



with k ^ {2n-l)T:/N and n = -A^/2 + 1, • • • , iV/2 where 
antiperiodic boundary conditions are chosen. The band- 
width \sW = 2nt. 



(iii) The two-dimensional Hubbard square lattice with 
nearest-neighbor hopping amplitude —t. The dispersion 
relation is given by 



e(k) ~ — 2i(cos k^ + cos ky) 



(9) 



with k = (2Tmx/L,2Tmy/L) where nx,ny = —L/2 + 
1, • • • , L/2 and N = L^. The bandwidth is W = 8t. 

The one-dimensional Hubbard model (i) is exactly 
solvable via the Bethe Ansatzll3|^d can easily be stud- 
ied using the real-space DMRGO Comparison with the 
exact solutions and the real-space method will pro- 
vide an opportunity to test the performance of the 
momentum-space DMRQ. The 1/r-Hubbard model (ii) 
is also exactly solvable,E3 but it is difficult to investi- 
gate with the real-space DMRG because the hopping 
is long-range and imaginary. For this model, the ad- 
vantage of the momentum-space approach is clear: one 
need only change the real, diagonal dispersion e{k). 
We will therefore use this model to investigate how a 
substantial change in the range of hopping affects the 
momentum-space algorithm. The Hubbard model on a 
two-dimensional square lattice (iii) will be used to inves- 
tigate and compare the effects of dimensionality on the 
momentum-space method and on the real-space method. 



III. DMRG IN MOMENTUM SPACE 

In principle, the usual DMRgII can be applied di- 
rectly to the momentum-space representation of the 
Hubbard model (^). In the momentum-space approach 
each Bloch function with momentum k and spin a cor- 
responds to a lattice site. To perform the calculations 
presented in this paper, we have adapted a program orig- 
inally written by WhiteB This program predates (and 
thus does not use) some recent developments which can 
greatly improve the performanc&.p£ the DMRG such as 
the wavefunction transformationutJ, the use of-eompos- 
ite operatorsjj and of non-abelian symmetries.ta Never- 
theless, this program is highly optimized and allows us 
to carry out DMRG calculations keeping up to to = 4000 
density-matrix eigenstates on a workstation with 1GB of 
memory. Below and in the next section we discuss fea- 
tures of the momentum space DMRG which differ from 
the usual real space DMRG. 

The summation in the second term of Eq. (^) runs over 
products of operators. A straightforward implemen- 
tation of the DMRG algorithm requires calculating and 
keeping track of 0{N^) matrices representing the differ- 
ent products of operators. This represents a significant 
increase compared to the real space-approach which re- 
quires only a constant number of operators for the one- 
dimensional and order 0{L ~ Vn) matrices for theiwo- 
dimensional Hubbard model (|l|) , respectively. XiangEl has 
shown that is it possible to define so-called composite op- 
erators and thus reduce the number of operators which 
need to be kept to 67V. In White's program, internal 
sums over blocks are carried out to reduce the number of 
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operators to 0{N'^) rather than 0{N), and an efficient 
representation of operators with small sparse matrices is 
used. 

We explicitly use the conservation of the particle num- 
ber iVc = ^ J of the z-componcnt of the total spin 
Sz — (1/2) ^ <7 ?T.fccr, and of the total momentum 



K = 



k.a 



k nkrr mod 27r . 



(10) 



Momentum symmetry reduces the size of the effective 
Hilbert space by about a factor N and allows us to de- 
compose the matrix representations of operators into sev- 
eral smaller matrices. Therefore, the dimension of the 
effective Hilbert space for a given number m of density 
matrix eigenstates is smaller in the momentum-space ap- 
proach than in the real space approach. This should be 
kept in mind when comparing results obtained with both 
approaches: In general, m can be made larger for a given 
amount of computational effort in momentum space. 

In a one-dimensional system in real space, there is a 
natural ordering of the lattice sites. . In two dimensions, 
there is some choice in the ocderingjiirllj but reasonable 
choices yield similar results.li3 In momentum space it is 
not a priori clear how one should arrange the sites in 
the lattice.H Thus, we have tested several possibilities in- 
cluding random ordering. We have found that the order 
of sites should be carefully chosen in the momentum- 
space approach - the rate of convergence and the accu- 
racy strongly depends on the site order. Fundamentally, 
it seems that Bloch states which are strongly scattered 
by the Hubbard term in (^ should be arranged to be as 
close together as possible. For the one-dimensional and 
two-dimensional Hubbard models with nearest-neighbor 
hopping we use an energetic ordering in which the sites 
are arranged according to |£(fe) — ep], where ep denotes 
the Fermi energy in the non-interacting case (U = 0). 
For the 1/r-Hubbard model, the Fermi energy has no 
particular relevance for low-energy scattering processes 
and we have found that ordering the sites according to 
e{k) works best. 

We use the finite-system DMRG algorithm and per- 
form several sweeps through the lattice until the ground- 
state energy converges as in the real-space approach. 
Wilson's numerical RG method is used instead of the 
infinite system algorithm to build up the lattice during 
the initial iteration. For the next iterations we apply 
the usual blocking scheme with a superblock made of 
two sites and two blocks with at most m states each. In 
Ref. [g] the superblock was built using two blocks and a 
single site. According to Xiang, this single-site approach 
is faster than the usual blocking scheme. As discussed 
in Ref. ^, however, the single-site blocking scheme is not 
a robust method unless several states are targeted. As 
we always target a single state in our calculations [the 
ground state for some quantum numbers {Nc, Sz, K)], 
we use the two-site blocking scheme. 

We have observed that the DMRG has difficulty find- 
ing the ground state when the interaction U is not weak 



or for some particular choices of the quantum number 
{Nc, Sz, K). The DMRG sometimes seems to converge 
first to a state other than the ground state and only con- 
verges to the true ground state after many sweeps or 
after the number m of density-matrix eigenstates is in- 
creased. This behavior is marked by a rapid drop in the 
energy after a relatively large number of sweeps or at a 
high value of m. A similar problem has been reported 
with the real-spaceJ3MRG applied to two-dimensional 
fermion systems O'Ej Therefore, one should not rely on 
DMRG results obtained for a fixed number m of density- 
matrix eigenstates kept or a fixed number of sweeps, but 
one should investigate the behavior of the DMRG as a 
function of m and of the number of sweeps. 



IV. CONVERGENCE AND ACCURACY OF 
MOMENTUM SPACE DMRG 

In this section we discuss the convergence and accuracy 
of the momentum-space DMRG method applied to the 
Hubbard model (||). We have also applied the real-space 
DMRG method to the real-space representation (|I]) of 
the Hubbard Hamiltonian with periodic boundary con- 
ditions. This allows us to make comparisons of both 
DMRG methods in order to illustrate both their differ- 
ences and their common features. Our real-space DMRG 
program uses more advanced techniques and is better 
optimized than our momentum-space DMRG program. 
Therefore, we have chosen to present no comparison of 
computer CPU time and memory usage in this paper as 
they would be meaningless. 

As a measure of the DMRG precision we use the error 
in the ground-state energy per site 



AE{m) 



EbMRG ("l) ~ -E'o: 
Wt 



(11) 



where inexact is the exact ground-state energy [for par- 
ticular quantum numbers {Nc, Sz, K)] and i?DMRG('7i) 
is the corresponding DMRG energy obtained with m 
density-matrix eigenstates kept. The exact results -Ecxact 
are calculated using the Bathe Ansatz for the one- 
dimensional Hubbard modeliJ and are computed numer- 
ically using exact diagonalizatkm techniques for the two- 
dimensional Hubbard model £j For the 1/r-Hubbard 
model they pare derived from a conjectured effective 
Hamiltonian.O Consistency of the DMRG energies with 
the spectrum obtained from the effective Hamiltonian in 
turn confirms the conjecture. 

Figures |l|(a) and (b) show the ground-state energy 
error AE of the one-dimensional Hubbard model as a 
function of the number of density-matrix eigenstates m 
for several values of U/t. These results have been ob- 
tained on 16-site lattices with periodic boundary condi- 
tions at half filling using the momentum-space [Fig. 0(a)] 
and the real-space [see Fig. ^(b)] approaches. The error 
AE' of the momentum-space DMRG method clearly in- 
creases with U/t. The error in the real-space DMRG 
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FIG. 1: DMRG error in the ground-state energy of the half-fiUed one-dimensional Hubbard model with 16 sites and periodic 
boundary conditions as a function of the number m of density-matrix eigenstates kept in (a) the momentum-space approach 
for U/t — 1 (circles), 2 (triangles), 3 (squares), 4 (reverse-triangles), and (b) the real-space approach for U/t = 1 (circles), 4 
(reverse-triangles), 8 (diamonds). 



increases with decreasing U/t for this half-filled sys- 
tem. In the momentum-space DMRG, the procedure 
becomes exact when the off-diagonal interaction terms 
vanish (i.e., at J7 = 0) and should be more accurate when 
they are small - it is a weak-coupling method. In con- 
trast, the real-space representation becomes exact (i.e., 
local) when i — > 0. It is important to note that this is not 
equivalent to the large U jt limit, in which the real-space 
DMRG does not become exact. The increase in accuracy 
with IJ jt shown here is specific to the half-filled insu- 
lator, in which the charge degrees of freedom become 
increasingly localized with increasing \J jt. In fact, in 
the one-dimensional system away from half band-filling, 
there is very little dependence ofihe convergence on IJ jt 
(for open boundary conditions). l3 In both approaches, 
Ai?(m) does not decrease exponentially as m increases, 
contrary to the behavior often reported for real-space 
DMRG calculations on one-dimensional systems with 
open boundary conditions. In Figs. |l|(a) and (b) we also 
see that the errors of the real-space approach are smaller 
and decrease more rapidly with m than the errors of the 
momentum-space approach. 

In Xiang's work,i the systematic convergence of the 
momentum-space method seems to break down when the 
interaction strength U approaches the bandwidth W . In 
particular, for the half-filled Hubbard model on a one- 
dimensional 16-site lattice with U jt = 4, the ground- 
state energy obtained from the DMRG, shown in Fig- 
ure ^, does not seem to converge smoothly toward the 
exact ground-state energy as the number m of retained 
density-matrix eigenstates is increased - the results are 
oscillatory and hard to extrapolate. While the origin of 
this irregular convergence in Xiang's data is unclear, one 
factor that is essential to consider is that the momentum 
of the ground state for a half-filled ring with 16 sites 
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FIG. 2: DMRG error in the ground-state energy as a function 
of the number m of density-matrix eigenstates kept for the 
16-site one-dimensional Hubbard model at half filling with 
U jt = 4, including our results for momentum K =^ (open 
circles) and A" = (filled circles) and Xiang's resuItsH for un- 
specified momentum (triangles). The dashed lines are guides 
to the eye. 



and periodic boundary conditions is X = tt because it 
is an open-shell configuration. Our results for K = ix 
do converge smoothly to the exact solution as a function 
of 1/m, as seen in Fig. ^ The lowest-energy state with 
K — 0, also shown in the figure, lies closer to Xiang's re- 
sults for the larger values of m, but converges smoothly 
to an energy that is clearly higher than the ground state. 
The deviation of Xiang's result could either be due to 
his ground state having a different momentum or due to 
convergence of the the DMRG to a state other than the 
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FIG. 3: DMRG errors in the ground-state energy at half fiUing as a function of the interaction strength U/W calculated with 
m = 2000 density-matrix eigenstates. The dashed lines are guides to the eye. (a) Real-space DMRG results for the Hubbard 
model on a one-dimensional lattice with 16 sites (circles) and 32 sites (squares), and on a two-dimensional 4x4 lattice (crosses) . 
(b) Momentum-space DMRG results for the one-dimensional Hubbard model with 16 sites (open circles) and 32 sites (open 
squares), the 1/r-Hubbard model with 16 sites (filled circles) and 32 sites (filled squares), and the two-dimensional Hubbard 
model on a 4 X 4 lattice (crosses). 



ground state, as discussed in the previous section. 



A. Dependence on model parameters 



For weaker interaction, U/t = 1, Xiang's results con- 
verge to the ground-state energy for K — ir, but lie 
significantly above our own results for the same num- 
ber of density-matrix eigenstates kept. For instance, 
the DMRG error in the ground-state energy AE{m) for 
m = 1200 is one order of magnitude smaller in our cal- 
culations than the value reported by Xiang. This also 
suggests an incomplete convergence of DMRG in Xiang's 
calculations even for weak coupling. However, it should 
also be kept in mind that in his work Xiang used a differ- 
ent superblock structure with a single site between two 
blocks. As a consequence, the dimension of the effective 
Hilbert space for a given number of states m is smaller 
in Xiang's calculations than in our calculations. This 
difference could be responsible for the better accuracy 
of our results in the weak coupling limit {U/t = 1) but 
cannot explain the discrepancy observed for intermediate 
couphng {U/t — 4). 

In our calculations we have often observed that the 
DMRG energy initially seems to converge towards a value 
larger than the exact ground-state energy. Upon fur- 
ther increasing the number of states m or the number 
of sweeps, it then converges to the exact result. In all 
the cases we present here, the momentum-space DMRG 
yields energies that do ultimately converge to the ex- 
act result, even for very large interaction strength U, 
although the rate of this convergence and the accuracy 
deteriorate rapidly as U increases. 



In this section, we discuss the dependence of the ac- 
curacy of the energy of momentum-space DMRG on the 
model parameters: the single-electron dispersion e(fc), 
the interaction strength U/W, the lattice dimensional- 
ity, and the band filling. Figures ^(a) and (b) show the 
ground-state energy error AE as a function of the in- 
teraction strength U /W for a fixed number m = 2000 
of density-matrix eigenstates. We again see that errors 
decrease in the real-space approach [Fig. ^(a)l but in- 
crease in the momentum-space approach [Fig. @(b)] for 
increasing U/W. In both cases errors increase with the 
system size and are larger in two dimensions than in one 
dimension for the same number of lattice sites. However, 
the dependence on system size and dimensionality of the 
momentum-space approach is clearly weaker than in the 
real-space approach. The lower precision of the DMRG 
in higher dimension is easily understood as a consequence 
of increasing off-diagonal coupling in the real-space ap- 
proach. A possible explanation for the slight decrease in 
accuracy with dimension in the momentum-space repre- 
sentation is that although the single-electron dispersion 
£{k) remains diagonal for any dimension, a larger propor- 
tion of single-electron states are close to the Fermi sur- 
face and are thus strongly scattered. In Figs. ||(a) and (b) 
one also sees that the precision of the real-space approach 
is generally better than that of the momentum-space ap- 
proach in the one-dimensional Hubbard model. The lat- 
ter becomes more accurate than the former for U/t 1 
only. In two dimension, however, the real-space approach 
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performs very poorly for periodic boundary conditions. 
The momentum-space approach yields better results for 
U < W = 8t. We finally note that DMRG errors seem 
to be affected only a small amount by the form of the 
dispersion s{k) in the momentum-space approach. In 
the real-space approach, however, changing the single- 
electron dispersion by introducing longer-range hopping 
lowers the DMRG performance very rapidly. In sum- 
mary, we find that the momentum-space approach is 
superior to the real-space approach for applications to 
translationally invariant systems with weak to interme- 
diate Coulomb interactions in two dimensions or on one- 
dimensional lattices with long-range hopping. 

Let us now consider the effects of the band filling. In 
Fig. ^ we show the error AE in the ground-state en- 
ergy as a function of band filling n = N^/N for the 
one-dimensional Hubbard model at C/ = At. In the 
momentum-space approach, the accuracy is worst at or 
near half filling and improves as the density decreases 
from n = 1. One cause of this effect is that the size 
of the Hilbert space is maximal at half filling and de- 
creases rapidly at large doping. This effect is magnified 
in the 16-site system relative to the 32-site system, as 
seen in Fig. ^, because a substantial proportion of the 
Hilbert space is retained in the diagonalization step at 
large doping. Another possible cause is the reduction in 
the effective strength of the electron-electron scattering 
as the system is doped away from half filling. The effec- 
tive interaction strength depends on the ratio of U and 
the density of states at the Fermi energy, which becomes 
smaller with doping, in weak coupling. As the effective 
interaction becomes smaller, the electrons become more 
localized in momentum space, which should be favorable 
for the convergence of the momentum-space DMRG. In 
the real-space approach. Fig. || shows that the error in 
the ground-state energy first increases as the system is 
doped slightly away from half filling, then decreases upon 
further doping. As discussed previously, the charge de- 
grees of freedom are localized in the half-filled insulator, 
leading to improved convergence for the real-space al- 
gorithm. For any finite doping, the system immediately 
becomes metallic, i.e., some charge degrees of freedom be- 
come delocalized, leading to a reduction in accuracy. As 
the system is doped further from half filling the reduction 
in the size of the Hilbert space leads to an improvement 
in accuracy, as in the momentum-space approach. 



B. Extrapolation to m ^ cxa 

DMRG calculations have a truncation error which is 
reduced by increasing the number m of density-matrix 
eigenstates kept.Cfd It is important to analyze the scal- 
ing of DMRG results as a function of m to estimate 
DMRG errors quantitatively. In real-space DMRG calcu- 
lations one generally observes that energy erroifS-|A£'(m) 
are proportional to the discarded weight Pm,lla'E2l pro- 
vided that the DMRG has converged to the right target 
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FIG. 4: DMRG error in the ground-state energy as function 
of band filling n in the one-dimensional Hubbard model for 
U/t — 4 on lattices with 16 sites (squares) and 32 sites (cir- 
cles). Open symbols represent the momentum-space DMRG 
results and filled symbols for the real-space DMRG results. 
Real-space DMRG errors for 16 sites are smaller than 10~^ 
and not shown. The number of density-matrix states is 
m — 1200 and m = 2000 for the 16-site and 32-site system, 
respectively. 



state. Here the discarded weight Pm is defined as the 
total weight (sum of the density-matrix eigenvalues) of 
the discarded density-matrix eigenstates, averaged over a 
sweep through all lattice sites in the finite-system DMRG 
algorithm. Thus, it is possible to extrapolate DMRG 
eigenenergies to the limit Pm — > 0. This procedure pro- 
duces extrapolated energies which are closer to the exact 
eigenenergies than the DMRG energies calculated for a 
given value of m. Moreover, the extrapolation yields re- 
liable quantitative error margins for the eigenenergies. 

In momentum space, however, we have found that the 
linear relationship between the energy errors AE{m) and 
Pm often does not hold, even for small P,„. In fact, 
we find that the dependence of P,„ on m can even be 
non-monotonic. An extrapolation to vanishing discarded 
weight P,n — > is therefore generally not possible. That 
such a non-monotonic behavior is found at all is surpris- 
ing at first glance because the discarded weight Pm of 
the exact density matrix for the system decreases mono- 
tonically with increasing m per definition. An exact den- 
sity matrix for the ground state of the Hubbard model 
can be calculated numerically in small systems using ex- 
act diagonalization. Using the results of such a calcu- 
lation on a iV = 12 lattice, we have found that the 
density-matrix eigenvalues Wi,i — 1,2, .. . appear to de- 
crease exponentially as a function of i in the asymptotic 
regime i ^ 1. As a consequence, the exact discarded 
weight Pm = J2iLi '^i ^^'^ the corresponding energy error 
AE{m) must also decrease exponentially with increasing 
TO. This is observed for density matrices calculated in 
the momentum-space approach as well as those obtained 
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in the real-space approach (both for open and periodic 
boundary conditions). Such an exponential falloff of the 
density-matrix jsigenvalues has been found for exactly 
solvable models.Ej In an actual DMRG calculation, how- 
ever, the density matrix is calculated self-consistently. 
Thus different density matrices can be obtained for dif- 
ferent m, and Pm can, in principle, be an arbitrary func- 
tion of TO, except for the condition limm^oo Pm = 0. We 
expect such effects to be largest where Pm is large and 
the self-consistently determined density matrix is a poor 
approximation to the exact one. 

In the momentum-space approach, the error in the 
DMRG energy AE{m) does not decrease exponentially 
with increasing m, but rather shows a power-law behav- 
ior in 1/to in the limit m ^ 1. We therefore extrap- 
olate the momentum-space DMRG results to vanishing 
truncation error by performing a least-squares fit of the 
DMRG energies i?DMRG (^) for several to to a n-th order 
polynomial in 1/to, 



Eat { - 

TO 



ai 



0-2 



Ec^ + — + — 



+ —. (12) 

TO" ^ ^ 



An extrapolated energy for vanishing truncation errors 
(to — > cx)) is given directly by the fit parameter Eao- The 
energy i?DMRG('7i) must be a monotonically decreasing 
function of to since the DMRG is a variational method 
and increasing m means increasing the variational sub- 
space dimension. Therefore, -E'dmrg(1/"t-) must sat- 



isfy the constraint 



dx 



> for X in the range 



< a; < 1/to' < 1, where m! is smaller than or equal to 
the smallest value of to used in the fit. An obvious con- 
sequence of this constraint is that the first-order term in 
the polynomial Ef^til/m) must satisfy ai > 0. We have 
found that the best fit under this constraint systemati- 
cally gives ai — 0. 

Fig. I shows the DMRG errors Ai? in the ground-state 
energies of the one-dimensional 32-site Hubbard model 
for U /t — 1 and 2 and the results of least-square fits 
to a fourth-order polynomial, n = 4 in Eq. (p^). The 
DMRG errors /^E{m) for the largest value of to used 
are 4.5 x 10"'' (m = 2000) and 5.3 x lO'^ (to = 2800) 
for U/t = 1 and U/t — 2, respectively. The accuracy is 
greatly improved by the polynomial fit and the to — > cxd 
extrapolation. The errors in the corresponding extrap- 
olated ground-state energies per site are 2.9 x 10~^ for 
U/t = 1 and 4.5 x 10"^ for U/t = 2. 

To illustrate the benefit of extrapolating DMRG ener- 
gies to vanishing truncation errors we now consider the 
quasi-particle gap Aqp of the one-dimensional Hubbard 



model at half filling. 

by 



The quasi-particle gap is defined 



= Eo{N + l:N) + Eo{N -1-N)-2Eq{N-N) , (13) 

where Eo{Nc]N) is the ground-state energy of a sys- 
tem with N sites and iVe electrons with minimal Sz, 
i.e., iV| = A^i or 7V| = A^X -I- 1. The charge con- 
jugation symmetry of the Hubbard model implies that 
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FIG. 5; DMRG error in the ground-state energy as a function 
of 1/m for the one-dimensional Hubbard model on a 32-site 
lattice. The lines are least-square fits to a fourth-order poly- 
nomial, n = 4 in Eq. (n2h . 



Eq{N -f 1; A) = E^{N - 1; A) + [/, so only one of these 
two energies need be computed. Some results for Aqp are 
shown in Table |. If one calculates the quasi-particle gap 
with the DMRG results for £;o(Ac; A^) obtained for fixed 
values of m, the magnitude of errors in Aqp fluctuates 
widely. The origin of this behavior has two competing 
sources. On the one hand, DMRG errors in the eigenen- 
ergies Eo{Ne;N) tend to be systematic for similar val- 
ues of A and A"o and cancel when calculating Eq. dlT 



Thus, the absolute error in Aqp can be smaller than 
the error in Eo{Nc; N) as seen in Table | for the case 
U = 2t. On the other hand, the eigenenergies are exten- 
sive quantities [i.e., Eo{Nc; N) scales with A^ for constant 
density N^/N], while Aqp is an intensive quantity (i.e., 
Aqp tends to a constant for increasing A and constant 
density Nc/N). Thus, even small but non-systematic 
errors in Eq{Nc', N) immediately result in much larger 
relative errors in Aqp. As a consequence, the values of 
Aqp (or similar physical quantities) obtained for a given 
number of density-matrix eigenstates kept might be ac- 
curate but there is considerable uncertainty about their 
accuracy. The extrapolation of the ground-state ener- 
gies Eo{Nc; A) to vanishing truncation errors allow us to 
eliminate this uncertainty and even to improve the pre- 
cision of our results. For instance, in Table | one can 
see that the errors in Aqp calculated with extrapolated 
ground-state energies are up to one order of magnitude 
smaller than for the largest value of to used (m = 2000). 



V. DISPERSION OF SPINON EXCITATIONS 

An advantage of the momentum-space method is that 
momentum-dependent quantities can be easily calcu- 
lated. In this section, we investigate the dispersion of the 
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FIG. 6: Spin excitations in the one-dimensional Hubbard model at half filling for U /t = 0.4. (a) Energy et{k;N) of the 
lowest triplet excitation calculated with DMRG as a function of the inverse system size 1/N for momenta fc = (filled circles), 
0.2n (filled squares), OAtt (filled triangles), O.Btt (open triangles), O.Stt (open squares), and n (open circles). Dashed lines are 
linear fits in 1/A'^. (b) Spinon dispersion Es{k) in the thermodynamic limit oo. The circles show the results obtained by 

extrapolation of the DMRG data in (a) [see Eq. ha] . The solid line is the exact result obtained in Ref. hi 



TABLE I: Ground-state energies and quasi-particle gap, Eq. (p_3), for the one-dimensional Hubbard model with = 32 sites. 



U/t 




m = 2000 


extrapolated 


exact 


1 


£;o(A^;A^)A 


-33.20078 


-33.21423 


-33.21515 




Eo{N+l-N)/t 


-32.64191 


-32.65757 


-32.65687 






0.12186 


0.11472 


0.11515 


2 


Eo{N; N)/t 


-26.80161 


-27.01970 


-27.01826 




EoiN + l-N)/t 


-25.61577 


-25.83275 


-25.83170 






0.37737 


0.37390 


0.37311 



spinon excitation Es{k) at half-band filling for both the 
one-dimensional Hubbard model with nearest-neighbor 
hopping and the l/r-Hubbard model; the spinon spec- 
trum is known exactly in both cases. 



The lowest spin-triplet excitation with momentum k 
in a system of size TV with a singlet ground state is given 

by 



et{k; N) = Eo{N,/2 + l,N,/2 -l,k + fco; N) - Eo{N,/2, N,/2, ko; N) , (14) 

where Eq(N^, N^, k; N) denotes the energy of the lowest state with (^i) up-spin (down-spin) electrons and total 
momentum k, and /cq is the momentum of the singlet ground state. In one-dimensional spin-1 /2 systems, a spin-triplet 
excitation is composed of two or more spin-1/2 spinous which are gapless in the thermodynamic limit. Therefore, the 
lowest-energy spin-one excitation allows us to map out the spinon dispersion as a function of momentum 

Es{k) ^ et{k; N ^ oo) ; < fc < tt . (15) 

In Fig. |6|(a), we show the lowest triplet excitation et{k; N) of the half-filled Hubbard model with nearest-neighbor 
hopping at U/t = 0.4 for several momenta fc and system sizes N up to 70 sites. The excitation energies scale 
approximately linearly with 1/N for all momenta but the slope varies considerably for the different fc. This finite-size 
scaling is readily understood from the exact result at [/ = for a closed-shell system, 

£t(fc;^) = 2isin(fc) + 8tsinQ')sin(^)sinQ-^') , (16) 
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where k = 27rn/N and n = 1, . . . , N/2. We therefore extrapolate the energies et{k; N) to the thermodynamic hmit 
using a hnear fit in 1/A^ to obtain the spinon dispersion EJk) though Eq. (|l^). The results for Es{k) are shown in 



Fig. ^(b) and compared with the exact spinon dispersionEa for U/t — 0.4. One sees that our numerical results for 
Es{k) agree very well with the analytical curve. In particular, et{k — 0,tt;N oo) — within the accuracy of the 
extrapolation. Note that because of the weak interaction U — OAt considered in this example, DMRG is very accurate 
and extrapolation of eigenenergies as a function of the number m of density-matrix states is not necessary. Thus, we 
have used fixed numbers of states m = 400, 800, 1200, 2000 for N = 10, 30, 50, 70, respectively. Actually, errors due to 
the infinite-system extrapolation are at least an order of magnitude larger than the DMRG error in the ground-state 
energy (per site), which is 2.7 x lO^"' in the worst case, = 70. _ 
In the 1/r-Hubbard model, the spinon spectrum at half band-filling is given byll3 




, AWUik-^) 2W{k-^) \ 
W^ + U^ ^ ^ + ^ ^ -U] , < yfc < TT 

TT TT I 



(17) 



in the thermodynamic limit. In order to form the low- 
est spin-triplet excitations, two spinon excitations are 
necessary. If they have momenta ki = k and k2 — 0, 
the spin-triplet excitation will have total momentum k 
and excitation energy et{k) = Es{k) with respect to the 
ground state, as i?s(fe = 0) = 0. In Fig. |j, this ana- 
lytical result for an infinite system is compared to our 
momentum-space DMRG results in finite systems of size 
TV = 24 and TV = 32 for U/t ^ 1. We observe an al- 
most perfect agreement, implying the absence of signifi- 
cant finite-size effects. This is a consequence of the linear 
dispersion and weak coupling considered here. It can be 
inferred from the results of Ref. ^ that the dispersion 
of spin-triplet excitations has no explicit dependence on 
the system size to first order in U/W 

et{k;N)^tk(^l-^^ . (18) 

Therefore, the finite-size corrections are of the order 
{1/N)(U/W)'^ <C 1 and are negligible in the results for 
U/t = 1 and > 24 presented in Fig. @. As in the 
nearest-neighbor hopping case, DMRG errors in the en- 
ergy are negligibly small compared to the spinon band- 
width because of the weak interaction used here. Thus 
in Fig. we show DMRG results for a fixed number 
m of density-matrix states (m = 800 for A^ = 24 and 
m = 1600 for A — 32) instead of results extrapolated to 
the m — > cx) limit. These results, along with the results 
for the one-dimensional Hubbard model discussed previ- 
ously, show that the low-lying spin-excitation spectrum 
can be accurately calculated using the momentum-space 
DMRG, at least in the weak-coupling limit. 

VI. MOMENTUM DISTRIBUTION 

Another quantity which is easily accessible to the 
momentum-space DMRG is the single-particle momen- 
tum distribution 




FIG. 7: Spinon spectrum Es{k) of the 1/r-Hubbard model 
for U/t = 1. DMRG results for finite systems with N ^ 24 
sites (squares) and A'' = 32 sites (circles). The solid curve is 
the exact result for an infinite system [Eq. (p^j. 



We have calculated n{k) in the ground state of the one- 
dimensional Hubbard model with nearest-neighbor hop- 
ping and of the 1 /r-Hubbard model at half band-filling 
using DMRG. Since DMRG truncation errors are typi- 
cally larger for quantities such as n(k) than for eigenen- 
ergies, it is crucial to examine the effect of varying m. 
The relative size of the finite-size effects is also impor- 
tant because one is interested in the behavior of n(k) in 
the thermodynamic limit. In the following, we compare 
our DMRG results to analytic results for n{k) on infinite 
lattices in the limit of both large and small interaction 
strength. 



A. Nearest— neighbor hopping 

In the ground state of the one-dimensional Hubbard 
model with nearest-neighbor hopping, the distribution 
n{k) is symmetric, n{—k) = n{k). Therefore, we show 
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results for < fc < tt only. 

We compare our DMRG results with an approxi- 
mate expression for n{k) proposed recently by Koch and 
Goedecker.LJ They make the ansatz that the real-space 
one-particle density matrix (cl^Cja) of the interacting 
system can be written as a product of the density matrix 
for the non-interacting system and an exponential decay 
factor that is a function of the particle-hole distance. 
The corresponding momentum distribution is 



1 1 / cos( 
"-KG(fc) = o + ~ arctan — - , 

2 TT \smh(7) / 



Here 1 denotes the decay length and is given by 



dx 



4{x)~Jl{x) 

1 + cxp(C/x/2i) 



(20) 



(21) 



where Jo (a:) and Ji{x) are Bessel functions of the first 
kind. In Ref. |23| it is found that exact diagonalization 
calculations on systems of up to 16 sites agree well with 
this form for U/t > 6. For smaller C//t, deviations are 
seen for wavevectors k ^ kp. 

Expansion of Eq. ( pT| ) in the strong-coupling limit 
yields 



In 



/2t7rln2 



and one recovers the perturbative resultc 



1 _ 21n(2)£(fc) 
''^''> - 2 U 



u:$>t. 



(22) 



(23) 



In Fig. g, we compare Eq. ( pOD for A; > with the 
momentum-space DMRG results for U/t — 10 and U/t — 
20 on small finite lattices. We use periodic boundary 
conditions {(p = 0) for systems with N = An + 2 sites 
and antiperiodic boundary conditions ((/> — n/N) for sys- 
tem with N = An sites (where n is an integer) to en- 
sure closed-shell configurations. The momentum-space 
DMRG calculations agree reasonably well with the ana- 
lytical result, but we note that the agreement becomes 
less good with increasing system size for U/t — 20. An 
analysis of the behavior for different m shows that our 
DMRG results underestimate n{k) for |fc| < fcp = tt/2 
and overestimate n(k) for |fc| > /cp- Therefore the devi- 
ations seen in Fig. are a consequence of DMRG errors 
which become larger for the larger system sizes and we 
would expect all results to lie on the analytic curve as in 
Ref. ^ in the m — ^ oo limit. _ . 
In the limit of weak coupling, 7 is given hyEB 



7 « 7C(3) 



U 
SttI 



(24) 



where ((z) is the Riemann Zcta-function [C(3) ~ 1.2]. 
In Fig. ^ we compare Eq. ( pO[ ) with our DMRG results 
for U/t = 1 and U/t = 2 on a 70-site lattice. We su- 
perimpose results with different phases < 4> < 2tt/N 




FIG. 8: Ground state momentum distribution function n{k) 
of the one-dimensional half-filled Hubbard model ai U/t — 20 
for system sizes A*' = 12 (crosses). A'' = 14 (triangles) and 
A'' = 16 (circles) keeping m = 2000 states, and at f//t = 10 
for TV = 16 sites (squares) keeping m = 1200 states. The lines 
correspond to the ansatz (12(1). 
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FIG. 9: Single-particle momentum distribution of the one- 
dimensional half-filled Hubbard model on a 70-site lattice for 
different phases 0. The circles are for U/t = 1 (with m = 800 
states kept) and the triangles for U/t = 2 (with m — 1200 
states kept). The lines correspond to Eq. ( |20[ ) for U/t = 1 
(solid) and U/t = 2 (dashed). The inset shows a blowup of 
the region [■k/2 - O.OIIStt, n/2 + O.OUStt] around the Fermi 
momentum. 



to obtain better resolution in the vicinity of the Fermi 
momentum kp — tt/2. Note that we can do this only in 
the insulating phase where the dependence of the ground 
state on the IjjjW-se (j) is negligible for small (j) because the 
Drude weightErEI (which is proportional to d'^Eo/d(f>^\^^o) 
vanishes. Unlike the limit of strong coupling, the agree- 
ment between our DMRG data and Eq. (^0|) is not good 
near fcp, as seen in the inset of Fig. |9| The first deriva- 
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tive of n{k) at A: = fcp is n'{kY) ~ —64 at UJt = 1 
and n'(fcF) « -10.3 at U/t ^ 2, whereas Eq. (|l|) yields 



'KG 



-23 at [//t = 1 and t^kgC^f 



-5.0 at _ 



U/t = 2. Although the ansatz (^0|) correctly describes 
the overall shape of n{k) deep in the insulating regime 
[U/t > 6), we find that it does not quantitatively recover 
the U/t and |A: — A:f| ^0 scaling limit. 

The behavior of the slope in the scaling limit can be 
understood by examining the Green function for the sine- 
Gordon model.tj The Fouriet, transform of the equal- 
time Green function becomescJ 



nitik > 0) 



arctan 



vi,ik~Tr/2) 

Aqp/2 



(25) 



where Uh = 2i is the holon velocity and Zq = 0.9219 has 
been calculated in Ref. The slope at A: = /cp = 7r/2 is 
then given by 

4Zo 



n'f^{k = kp) 



TrAn 



which has the values n'f^{k = kp) 
and fT.ft(fc — k-p) = —6.8 at U/t 



= 2 



(26) 

233 at U/t = 1 
when, the exact 



values for Aqp in the thermodynamic limitE3 are used. 
In Fig. |l^ we show the scaling with inverse system size 
of the DMRG results for l/\n'{k = kp)]. For U/t = 1, 
a 1/iV — > extrapolation (using a linear form) yields 
n'{k — kp) ~ —207, which agrees to within the fit errors 
with the field-theoretic result, corroborating both the 
field-theoretic and DMRG results in this regime. For 
U/t = 2, we obtain n'{k = kp) = —12, whose absolute 
value is substantially larger (i.e., outside our estimate 
for the error) than that of the field-theoretic result. We 
believe that this is because U/t = 2 is large enough so 
that the field-theoretic approximation to the momentum 
distribution of the Hubbard model begins to show sig- 
nificant deviations from the lattice result. This is sup- 
ported by the fact that the field-theoretic formula (^) 
correctly reproduces the strong-coupling form (p3|), but 
with a prefactor 4Zo/7r « 1.2 rather than 41n2 « 2.8, in- 
dicating that the field-theoretic formula underestimates 
the slope at larger values of U. 



B. 1/r— hopping 

In the ground state of the 1 /r-Hubbard model at half 
filling, n{—k) = 1 — n{k). The half-filled model describes 
a metal for [/ < W — 2Tit and an insulator for U > W. 

The momentum distribution in the 1 / r-Hubbard mod- 
el can be calculated in perturbation theory. Fip, large 
coupling U ^ t, we use Takahashi's approached. To 
this end, we need the spin-spin correlation function in 
the Gutzwiller-projected paramagnetic Fermi-sea at half 
band-filling, whichjs the ground state of the Haldane- 
Shastry spin chain£3 The spin-spin correlatipii function 
for the Haldane-Shastry model is given byHEJ 

1^(5..5,,„> = (-1)^^ , (.^0), (27) 
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FIG. 10: Inverse of the slope at the Fermi wavevector 
l/\n'{k — tt/2)\ as a function of inverse system size 1/A'' for 
U/t — 1 (open circles) and U/t = 2 (open triangles). The 
lines are linear fits to the DMRG data. The filled symbols 
are the corresponding field-theoretic results and__the crosses 
are calculated from the Koch-Goedecker ansatz 



where Si(a;) denotes the sine-integral. The momentum 
distribution for large U/t is then 



n{k) 



1 1 



, sin(fcr) / 3Si(7rr) 



7rr 



1-- ( l-31n 



(28) 



Note that although the momentum distribution is con- 
tinuous at kp = 0, we can expect to observe a large ap- 
parent jump in numerical simulations of finite systems 
even at large U/t because of the sizable logarithmic term 
in Eq. 

For small couplings, standard perturbation theory 
in U/t is applicable because the model reduces to a pure 
54 model in the conformal limit .E£l It turns out that the 
Gutzwiller wave function becomes exact in the small- 
coupling limit so that the momentum distribution for 
U <^t becomesE^ 



n{k) 



1 - (C/VVF2)/(fc) for -TT < fc < 
{U'^/W'^)J{k) for < A: < TT ^^^> 



with f{k\^ 3/16 - (1/4 - |fc|/(27r))2 
Figure 



11(a) displays the momentum distribution func- 
tion for the half-filled 1 /r-Hubbard model for [//t = 20 
on an A'^ = 12 lattice compared with the strong-coupling 
result (|28|), and for U /t — A on an N = 16 lattice com- 
pared with the the weak-coupling result (29). The agree- 
ment is very good because the DMRG errors are negli- 
gible for such small systems and because the finite-size 
effects are small, at least on the scale of the figure. 

The finite-size effects in weak coupling are more visible 
m Fig. |n|(b). One can see that the perturbation theory, 
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FIG. 11: Momentum distribution function for the 1/r-Hubbard model at half filling with m = 2000 states kept, (a) DMRG 
results on an = 16 lattice for U/t = 4 and an A'' = 12 lattice for U/t = 20. The solid and dashed lines are from the 
perturbative first-order result in t/U, Eq. (^8|), and the second-order result in U/t, Eq. (|29|), respectively, (b) The weak- 
coupling regime for U /t = 1. 2 . 3 (from bottom to top). The crosses, open triangles and filled circles denote A*' = 16, 24, 32. The 
solid line results from Eq. (|29[). 



Eq. (^), agrees well with the numerical results at all sys- 
tem sizes for U/t = 1. However, deviations that become 
larger with system size can be seen for U/t = 2. The re- 
sults for the smaller systems incorrectly suggest that Eq. 
( p9|) applies perfectly and that finite-size effects are ab- 
sent. This effect is even more marked at U/t = 3 where 
the analytic weak-?7 result agrees almost perfectly with 
the numerical data for N = 16, but the data for N = 24 
reveal that the agreement would be worse in the ther- 
modynamic limit. Apparently, the finite-size effects are 
approximately the same order and sign as the higher- 
order corrections in U /t here. 

VII. DISCUSSION AND OUTLOOK 

In this work, we have examined in detail the applica- 
tion of the Density Matrix Renormalization Group to the 
momentum-space representation of the Hubbard model. 
We have treated three different dispersion relations cor- 
responding to the one-dimensional chain with nearest- 
neighbor hopping, the one-dimensional chain with hop- 
ping that decays as 1/r, and the two-dimensional square 
lattice with isotropic nearest-neighbor hopping. While 
the one- and two-dimensional nearest-fiieighbor hopping 
cases were treated previously by Xiangj3 here we have ex- 
tended the scope of his results and have addressed some 
issues raised by his work. In particular, we hajzE taken 
advantage of the Bethe Ansatz exact solutionjl3 which 
yields the exact ground-state energy to within machine 
precision for all system sizes, to make more extensive 
studies of the convergence for the one-dimensional model 
with the number of density-matrix eigenstates kept, m. 
We have examined the effect of system size, interaction 



strength and band filling at up to m = 4000. 

For all parameters and models, we have found system- 
atic variational convergence with m to the true ground 
state; this does not seem to clearly occur in Xiang's, at 
least in his U/t — A results for the one-dimensional sys- 
tem. However, the convergence seems to be slower than 
exponential for the range of m accessible to us. The 
accuracy decreases regularly with system size when the 
other parameters are fixed, as one finds for the real-space 
DMRG. The accuracy also becomes rapidly worse with 
interaction U, as also found by Xiang. Our more exten- 
sive set of interaction strengths indicates that the behav- 
ior is quite regular as a function of U ; convergence does 
not break down at a particular U value. 

In the one-dimensional Hubbard model with periodic 
boundary conditions, a comparison with the real-space 
DMRG applied to the same system with the same m 
indicates that the accuracy of the real-space method is 
better for U/t > 1. The dependence of the accuracy on 
band filling is weak in momentum space, except on small 
system sizes for which the proportion of the Hilbert space 
kept changes drastically with filling. 

For the two-dimensional system, we have restricted 
ourselves to the half-filled 4x4 system, the largest 
for which exact diagonalization data are available, and 
have compared with the real-space DMRG. We find that 
the momentum-space method is more accurate than the 
real-space method (for the same m) when t/ < 8t = W. 

One crucial issue raised by Xiang is the dependence of 
the accuracy on dimensionality or, relatedly, the range 
of the hopping. He speculated that such effects would 
be smaller for the momentum-space DMRG than for the 
real-space DMRG, a speculation that we have confirmed 
here. We emphasize however that the choice of values 
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of interaction strength which are compared is important. 
It is our opinion that a reasonable choice is the interac- 
tion divided by the bandwidth, U/W (at identical filling 
and number of lattice sites). The bandwidth sets the 
energy scale for many physical phenomena and also is 
relevant to the strength of the coupling in perturbation 
theory. (A possibly useful alternative might be the ra- 
tio of the interaction strength and the density of states 
at the Fermi energy, the coupling parameter in weak- 
couphng perturbation theory.) For given U/W (as op- 
posed to given U/t), the accuracy of the momentum- 
space DMRG is lower in two dimensions than in one, an 
effect which becomes smaller as the interaction becomes 
larger. In one dimension, changing to the longer-range 
1 /r hopping has little effect on the accuracy at weak U/t, 
although the accuracy does become somewhat worse as 
U/t is increased. We therefore conclude that while the 
performance of the momentum-space DMRG is less de- 
pendent on range of the hopping or dimensionality than 
the real-space DMRG, there is still some effect. 

While the ground-state energy is an important in- 
dicator of convergence, it is not a particularly useful 
quantity in determining the physical behavior of a sys- 
tem. We have therefore examined a number of other 
quantities that are easily accessible to the momentum- 
space DMRG, which also yield useful physical informa- 
tion. Gaps formed from differences in energies pro- 
vide important information about the excitation spec- 
trum. We have examined the single-particle gap as well 
as the momentum-dependent triplet gap for both one- 
dimensional models at half filling. For the single-particle 
gap, extrapolation in 1/m is crucial in obtaining consis- 
tent accuracy because of cancellation of variational errors 
and because the gap is an intensive quantity obtained by 
subtracting extensive energies. We have found that a 
direct extrapolation using a polynomial in 1/m is the 
best method because the error in the ground-state en- 
ergy is not well correlated with the weight of the dis- 
carded density-matrix eigenvalues, unlike in the real- 
space DMRG. For the spin excitation spectrum, we have 
treated parameter values for which extrapolation in m is 
unnecessary and found that the finite-size effects are sub- 
stantial for the nearest-neighbor-hopping chain. For the 
l/i — Hubbard model, finite-size effects were quite small. 
In both cases, the size-extrapolated spectrum agrees well 
with exact results. 

Finally, we have examined the momentum distribution 
function. For the nearest-neighbor chain, we corpeare 
with an analytical ansatz of Koch and Goedecker.H At 
strong coupling, we find very good agreement aside from 
deviations due inaccuracy of the DMRG results. At weak 
interaction, the DMRG results agree with field-theoretic 
resultsJ2j whereas there is significant deviation from the 



Koch-Goedecker ansatz. For the l/r-Hubbard model, 
agreement with weak and strong coupling results is good, 
although finite-size corrections with the same sign as 
higher order terms in U/t provide better agreement for 
small system sizes than is justified in the thermodynamic 
limit. 

In summary, the momentum-space DMRG can be a 
useful tool for the Hubbard model at weak to interme- 
diate coupling. While it is competitive with the real- 
space DMRG only at quite weak coupling for the one- 
dimensional model (even with periodic boundary condi- 
tions), it is competitive up to significantly stronger cou- 
pling for longer-range hopping or in two dimensions. It 
should be noted that the accuracy of the momentum- 
space DMRG is generally significantly lower than that 
of the more favorable cases for the real-space DMRG, 
and can be considered a "numerically exact" method only 
with reservations. It can, however, be a useful variational 
method where no more exact methods are available if its 
limitations are well understood. The ease of calculation 
of momentum-dependent quantities is very useful - such 
quantities are often not available on large systems even 
for well-understood models. For these quantities, care 
must be taken with respect to accuracy and finite-size ef- 
fects, but we have found them to be well-behaved within 
these limitations. 

The momentum-space DMRG code used here is far 
from maximally optimized. With similar optimizations 
as used in the real-space program such as the use of 
wavefunction transformations to improve theiiiitial guess 
for the target state in the diagonalization,tj it should 
be possible to keep significantly more density-matrix 
eigenstates for given numerical efi'ort in the momentum- 
space method than in the real-space method applied to 
the same system. Use of such a better optimized pro- 
gram might increase the range of applicability of the 
momentum-space DMRG somewhat. In addition, since 
the Hubbard model has a local interaction in real-space, 
the interaction is quite non-local in momentum space. 
The momentum-space DMRG could quite possibly be 
better suited to models with longer range interaction in 
real space, corresponding to more local interactions in 
momentum space. Such directions would certainly be 
worth exploring in future work. 
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